/*
  File:    Pressure.cpp
  Purpose: Solves the pressure equation. 

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>

using namespace std;

#include "../ImprovedMath.h"
#include "../Matrix.h"
#include "GhostFluidMethod.h"
#include "Curvature.h"
#include "GaussianFilter.h"
#include "DiffusionSolverSteadyState.h"
#include "AdvectTumor.h"
#include "NecrosisLevelSets.h"
#include "../LogfileEntry.h"
#include "Oxygen.h"

#include "../AngiogenesisInterface.h"

#include "../ParseText.h"

#include "Pressure.h"


char PressureBoundaryConditionType = 'D';



double MobilityMinECM_Cutoff;
double MobilityMaxECM_Cutoff;

Polynomial MobilityECM_Polynomial;
bool MobilityECM_PolynomialSet = false;

bool ShiftPressureToMatchAverageAfterUpdate = false;
bool PressureAverageIsSet = false;
double PressureAverage = 0.0;

double TumorCenterX;
double TumorCenterY; 

double mu_out;
double mu_in;

double o2_ext;
double o2_H;
double o2_N;

double R_far;
double R;
double R_H;
double R_N;

double ExternalRadialPressure_const = 0.0;
double ExternalRadialPressure_logterm = 0.0;

bool UseFarFieldRadialDirichletCondition = false;
double FarFieldRadialPressure_DissipationDistance = 100.0;

bool ShiftPressureToMatchAverage( void )
{
 // exit if the option is not set;
 
 if( !ShiftPressureToMatchAverageAfterUpdate )
 { return false; }
 
 extern Matrix Pressure;
 
 // calculate the current average pressure
 
 double CurrentPressureAverage = 0.0;
 for( int i=0 ; i < Pressure.TellRows() ; i++ )
 {
  for( int j=0 ; j < Pressure.TellCols() ; j++ )
  {
   CurrentPressureAverage += *Pressure(i,j);
  }
 }
 CurrentPressureAverage /= ( Pressure.TellRows() * Pressure.TellCols() * 1.0 );

 // Check to see if the (initial) pressure average has ever been 
 // set. If not, then overwrite it with the current pressure average. 
 
 if( !PressureAverageIsSet )
 {
  PressureAverage = CurrentPressureAverage;
  PressureAverageIsSet = true;
 }
 
 // shift the pressure;
 
 double PressureShift = PressureAverage - CurrentPressureAverage;
 
 cout << endl << "Shifting pressure from current average " << CurrentPressureAverage 
      << " to new average " << PressureAverage << " ... " << endl << endl;
 
 for( int i=0 ; i < Pressure.TellRows() ; i++ )
 {
  for( int j=0 ; j < Pressure.TellCols() ; j++ )
  {
   *Pressure(i,j) += PressureShift;
  }
 }
 
 return true;
}

bool CreateMobilityECM_Polynomial( void )
{
 if( MobilityECM_PolynomialSet )
 { return false; }

 extern char* ParameterFile; 
 
 char* szTemp;
 szTemp = new char [100];
 RetrieveStringParameter( ParameterFile, "PressureBoundaryConditionType", szTemp );
 PressureBoundaryConditionType = szTemp[0];
 if( PressureBoundaryConditionType != 'D' && PressureBoundaryConditionType != 'd' &&
     PressureBoundaryConditionType != 'N' && PressureBoundaryConditionType != 'n' && 
	 PressureBoundaryConditionType != 'E' && PressureBoundaryConditionType != 'e' )
 { PressureBoundaryConditionType = 'D'; }
 
 ShiftPressureToMatchAverageAfterUpdate = (bool) RetrieveDoubleParameter(ParameterFile,"ShiftPressureToMatchAverageAfterUpdate");
 UseFarFieldRadialDirichletCondition = (bool) RetrieveDoubleParameter(ParameterFile,"UseFarFieldRadialDirichletCondition");
 FarFieldRadialPressure_DissipationDistance = RetrieveDoubleParameter(ParameterFile,"FarFieldRadialPressure_DissipationDistance");
 if( FarFieldRadialPressure_DissipationDistance < 1 )
 { FarFieldRadialPressure_DissipationDistance = 100.0; }
 
 extern double MinimumMobility;
 MinimumMobility = RetrieveDoubleParameter(ParameterFile,"MinimumMobility");
 if( MinimumMobility < 0.001 )
 { MinimumMobility = 0.001; }

 extern double MaximumMobility; 
 MaximumMobility = RetrieveDoubleParameter(ParameterFile,"MaximumMobility");
 
// extern double MobilityMinECM_Cutoff;
 MobilityMinECM_Cutoff = RetrieveDoubleParameter(ParameterFile,"MobilityMinECM_Cutoff");

// extern double MobilityMaxECM_Cutoff;
 MobilityMaxECM_Cutoff = RetrieveDoubleParameter(ParameterFile,"MobilityMaxECM_Cutoff");
 if( MobilityMaxECM_Cutoff < 1e-8 ) 
 { MobilityMaxECM_Cutoff = 1.0; } 

 MobilityECM_Polynomial.SetDegree(3);
 MobilityECM_Polynomial.SetVariable('E');
 
 for( int i=0 ; i <= MobilityECM_Polynomial.GetDegree() ; i++ )
 { MobilityECM_Polynomial.SetCoefficient( i , 1.0 ); }

 Polynomial Prime = Derivative( MobilityECM_Polynomial );
 
 Matrix Coeff = zeros( MobilityECM_Polynomial.GetDegree() + 1 , 
                       MobilityECM_Polynomial.GetDegree() + 1 );
 Matrix RHS = zeros( 1, MobilityECM_Polynomial.GetDegree() + 1 );

 for( int j=0; j <= MobilityECM_Polynomial.GetDegree() ; j++ )
 { *Coeff(0,j) = MobilityECM_Polynomial.GetCoefficient(j)*IntPow( MobilityMinECM_Cutoff , j ); }
 *RHS(0) = MaximumMobility;
 
 *Coeff(1,0) = 0.0;
 for( int j=1; j <= MobilityECM_Polynomial.GetDegree() ; j++ )
 { *Coeff(1,j) = Prime.GetCoefficient(j-1)*IntPow( MobilityMinECM_Cutoff , j-1 ); }
 *RHS(1) = 0.0;

 for( int j=0; j <= MobilityECM_Polynomial.GetDegree() ; j++ )
 { *Coeff(2,j) = MobilityECM_Polynomial.GetCoefficient(j)*IntPow( MobilityMaxECM_Cutoff , j ); }
 *RHS(2) = MinimumMobility;
 
 *Coeff(3,0) = 0.0;
 for( int j=1; j <= MobilityECM_Polynomial.GetDegree() ; j++ )
 { *Coeff(3,j) = Prime.GetCoefficient(j-1)*IntPow( MobilityMaxECM_Cutoff , j-1 ); }
 *RHS(3) = 0.0; 

 Matrix FinalCoefficients = GaussianEliminationPartialPivoting( Coeff, RHS  );
 
 for( int i=0; i <= MobilityECM_Polynomial.GetDegree() ; i++ )
 { MobilityECM_Polynomial.SetCoefficient( i , *FinalCoefficients(i) ); }
 Prime = Derivative( MobilityECM_Polynomial );

 cout << "Mobility ECM polynomial details: " << endl 
      << "-------------------------------" << endl
	  << "Mu(ECM) = "; 
 MobilityECM_Polynomial.display();
 cout << endl;
// Prime.display();

 
 cout << "At " << MobilityMinECM_Cutoff << " : Mu(ECM) = " 
      << MobilityECM_Polynomial(MobilityMinECM_Cutoff) << " Mu'(ECM) = " 
	  << Prime(MobilityMinECM_Cutoff) << endl;
 
 cout << "At " << MobilityMaxECM_Cutoff << " : Mu(ECM) = " 
      << MobilityECM_Polynomial(MobilityMaxECM_Cutoff) << " Mu'(ECM) = " 
	  << Prime(MobilityMaxECM_Cutoff) << endl << endl;
 
 MobilityECM_PolynomialSet = true;

 return true;
}


double MobilityFunction(double ecm, double o2, double S)
{
 extern bool FunctionalMobility; 
 if( !FunctionalMobility )
 {
  if( S > 0.5 )
  { return 0.001; }
  return 1.0; 
 }
 
 CreateMobilityECM_Polynomial();
 
 extern double MobilityMinECM_Cutoff;
 extern double MobilityMaxECM_Cutoff;
 
 if( ecm < MobilityMinECM_Cutoff )
 { ecm = MobilityMinECM_Cutoff; }
 if( ecm > MobilityMaxECM_Cutoff )
 { ecm = MobilityMaxECM_Cutoff; }
 
 return MobilityECM_Polynomial(ecm)*(1.0-S) + 0.001*S; 
 
 extern double MinimumMobility;
 extern double MaximumMobility;

 
 
 extern double MobilityECM_ParameterA;
 extern double MobilityECM_ParameterB;
 
 
 double output = 1.0/( MobilityECM_ParameterA*ecm + MobilityECM_ParameterB + S*1000);
 if( output > MaximumMobility && S < 0.5 )
 { return MaximumMobility; }
 if( output < MinimumMobility && S < 0.5 )
 { return MinimumMobility; }
 return output;
}

bool UpdateMobility( void )
{
 extern double _EPS;
 
 extern int Xnodes;
 extern int Ynodes;

 extern Matrix ECM;
 extern Matrix Structure;
 extern Matrix Mobility;
 extern Matrix Oxygen;

 for( int i=0 ; i < Xnodes ; i++ )
 {
  for( int j=0; j < Ynodes ; j++ )
  {
     *Mobility(i,j) = MobilityFunction(*ECM(i,j),*Oxygen(i,j),*Structure(i,j));
	 //if(i==Xnodes-1&&j==Ynodes-1)
	 //{
	  // cout<<"here"<<i<<" "<<j;
      //getchar();
    // }

  }
 }

 extern Matrix X;
 extern Matrix Y;
 extern double dx;
 
 extern double MinimumMobility;
 
 extern double FunctionSmoothingLevel;
 Mobility = GaussianFilter(X,Y,Mobility,FunctionSmoothingLevel*dx);
 
 for( int i=0; i < Xnodes ; i++ )
 {
  for( int j=0; j < Ynodes ; j++ )
  {
   if( *Mobility(i,j) < _EPS )
   { 
    *Mobility(i,j) = 0.001;
   }
  } 
 }

 
 return true;
}

double PressureBC( double x )
{ return 0.0; }

double PressureJumpFunction( double x, double y )
{
 extern Matrix X;
 extern Matrix Y;
 extern Matrix Phi;
 extern double dx;
 extern double TumorSurfaceTensionConstant;
 extern double TumorAggressivenessG;
 
 GeometricDataPacket GDP = AdaptiveGeometricData(X,Y,Phi,x,y);
 return (TumorSurfaceTensionConstant/TumorAggressivenessG)*GDP.Curvature;
}

double PressureNormalJumpFunction( double x, double y )
{
 extern Matrix X,Y;
 extern double dx,dy;
 int i = (int) floor( (x-*X(0))/dx );
 int j = (int) floor( (y-*Y(0))/dy );
 
 double ThetaX = (x-*X(i))/dx;
 double ThetaY = (y-*Y(j))/dy;
 
 extern Matrix TaxisJump;
 
 // use a linear interpolation 
 
 if( ThetaX > 0.01 )
 { return *TaxisJump(i,j) * (1.0-ThetaX) + *TaxisJump(i+1,j) * ThetaX; }
 else
 { return *TaxisJump(i,j) * (1.0-ThetaY) + *TaxisJump(i,j+1) * ThetaY; }

 return 0.0; 
}

bool UpdatePressure( void )
{
 extern double _EPS;
 
 extern int Xnodes;
 extern int Ynodes;
 extern Matrix X;
 extern Matrix Y;
 extern double dx, dy;
 
 extern Matrix Phi;
 extern Matrix PhiNecrotic;
 extern Matrix PhiQuiescent;
 
 extern Matrix Oxygen;
 extern Matrix Mobility;
 extern Matrix Pressure;
 
 extern double ApoptosisRate;
 extern double TumorNecrosisRateGN;
 extern double TumorProliferationConstant;

 Matrix PressureSource = zeros(Xnodes,Ynodes);
 Matrix PressureSourceU = zeros(Xnodes,Ynodes);
 
 extern Matrix ECMtaxisCoefficient, ECM;
 
 UpdateECMtaxisCoefficient();
	extern Matrix TaxisJump;
 TaxisJump = zeros(Xnodes,Ynodes);
 
 extern double FunctionSmoothingLevel;
 Matrix ECMsmoothed = GaussianFilter(X,Y,ECM,FunctionSmoothingLevel*dx); 
 
 for( int i=1 ; i < Xnodes-1 ; i++ )
 {
  for( int j=1 ; j < Ynodes-1 ; j++ )
  {
   bool Done = false;
   if( *PhiNecrotic(i,j) <= _EPS )
   {
    *PressureSource(i,j) = -TumorNecrosisRateGN; 
    Done = true; 
   }
   if( *PhiQuiescent(i,j) <= _EPS && !Done )
   { Done = true; }
   if( *Phi(i,j) <= _EPS && !Done )
   {
    *PressureSource(i,j) =  TumorProliferationConstant*( *Oxygen(i,j) - ApoptosisRate);
   }
   
   double ECMtaxisXr   = 0.5*( *ECMtaxisCoefficient(i+1,j) + *ECMtaxisCoefficient(i,j) );
   double ECMtaxisXell = 0.5*( *ECMtaxisCoefficient(i-1,j) + *ECMtaxisCoefficient(i,j) );
   double ECMtaxisYr   = 0.5*( *ECMtaxisCoefficient(i,j+1) + *ECMtaxisCoefficient(i,j) );
   double ECMtaxisYell = 0.5*( *ECMtaxisCoefficient(i,j-1) + *ECMtaxisCoefficient(i,j) );
   
   double ECM_taxis_contribution = ( ECMtaxisXr*(*ECMsmoothed(i+1,j)-*ECMsmoothed(i,j)) -
                                     ECMtaxisXell*(*ECMsmoothed(i,j)-*ECMsmoothed(i-1,j))  )/square(dx) 
                                 + ( ECMtaxisYr*(*ECMsmoothed(i,j+1)-*ECMsmoothed(i,j)) - 
								     ECMtaxisYell*(*ECMsmoothed(i,j)-*ECMsmoothed(i,j-1))  )/square(dy);
 								 
   *PressureSourceU(i,j) -= ECM_taxis_contribution; 
   
   GeometricDataPacket GDP = AdaptiveGeometricData(X,Y,Phi,*X(i),*Y(j));

   *TaxisJump(i,j) = *ECMtaxisCoefficient(i,j)* (  GDP.NormalX*((*ECMsmoothed(i+1,j)-*ECMsmoothed(i-1,j))/(2.0*dx))
											      +GDP.NormalY*((*ECMsmoothed(i,j+1)-*ECMsmoothed(i,j-1))/(2.0*dy)) );
  }
 }
 
 UpdateMobility();

 int MaxIterations = 20000;
 double Tolerance = IntPow(dx,4);
 char bctype = 'D'; // 'N';
 
 //bctype = PressureBoundaryConditionType;

 double (*PressureBC_left)(double) = PressureBC;
 double (*PressureBC_right)(double) = PressureBC;
 double (*PressureBC_bottom)(double) = PressureBC;
 double (*PressureBC_top)(double) = PressureBC;

 
 if( UseFarFieldRadialDirichletCondition )
 {
  cout << "Using radial far-field tumor pressure for Dirichlet boundary condition ... " << endl;
  CalculateEquivalentCircularTumor();
  bctype = 'D';
  
  PressureBC_left = FarFieldRadialSolution_left;
  PressureBC_right = FarFieldRadialSolution_right;
  PressureBC_bottom = FarFieldRadialSolution_bottom;
  PressureBC_top = FarFieldRadialSolution_top;
 }
		  
 Pressure = GhostFluidMethod(X,Y,Phi,
			  Mobility, PressureSource,PressureSourceU, 
              PressureBC_left,bctype, PressureBC_right,bctype, 
              PressureBC_bottom,bctype, PressureBC_top,bctype, 
			  PressureJumpFunction, PressureNormalJumpFunction, 
			  Pressure, MaxIterations, Tolerance );
			  
 ShiftPressureToMatchAverage();	

/******

	 double v = (pow(5, 0.5) - 0.1)/(pow(5,0.5) - 1);
 
 
  for(int i=0; i< Xnodes; i++)
 {
 for(int j=0; j< Ynodes; j++)
 {
 double r = pow( pow((j-120)*1.0/200,2) + pow((i-100)*1.0/200, 2) , 0.5);
 if(r<= 0.1)
 {
 *Pressure(i, j) = 1;
 }
 else
 {
 *Pressure(i, j) =(v - r)*(v - r)*1.0/pow((v-0.1), 2);

 }
 // *TAF(i, j) = exp(-pow((j)*1.0/200,2)/0.45); 
 //cout<<*TAF(i, j) <<endl;

 }

 
 }





*****************/ 
 return true;
}

bool QuicklyUpdatePressure( void )
{
 extern double _EPS;
 
 extern int Xnodes;
 extern int Ynodes;
 extern Matrix X;
 extern Matrix Y;
 extern double dx;
 
 extern Matrix Phi;
 extern Matrix PhiNecrotic;
 extern Matrix PhiQuiescent;
 
 extern Matrix Oxygen;
 extern Matrix Mobility;
 extern Matrix Pressure;
 
 extern double ApoptosisRate;
 extern double TumorNecrosisRateGN;
 extern double TumorProliferationConstant;

 Matrix PressureSource = zeros(Xnodes,Ynodes);
 Matrix PressureSourceU = zeros(Xnodes,Ynodes);

 extern double QuiescenceOxygenLevel;
 extern double NecrosisOxygenLevel;
 
 for( int i=0 ; i < Xnodes ; i++ )
 {
  for( int j=0 ; j < Ynodes ; j++ )
  {
   bool Done = false;
   if( *PhiNecrotic(i,j) <= _EPS || *Oxygen(i,j) <= NecrosisOxygenLevel )
   {
    *PressureSource(i,j) = -TumorNecrosisRateGN; 
    Done = true; 
   }
   if( *PhiQuiescent(i,j) <= _EPS || *Oxygen(i,j) <= QuiescenceOxygenLevel )
   {
    Done = true; 
   }
   if( *Phi(i,j) <= _EPS && !Done )
   {
    *PressureSource(i,j) =  TumorProliferationConstant*( *Oxygen(i,j) - ApoptosisRate);
   }
  }
 } 
 
 UpdateMobility();

 int MaxIterations = 2000;
 double Tolerance = IntPow(dx,4);
 char bctype = 'D'; // 'N';
			  
 Pressure = GhostFluidMethod(X,Y,Phi,
			  Mobility, PressureSource,PressureSourceU, 
              PressureBC,bctype, PressureBC,bctype, 
              PressureBC,bctype, PressureBC,bctype, 
			  PressureJumpFunction, PressureNormalJumpFunction, 
			  Pressure, MaxIterations, Tolerance );
			  
 ShiftPressureToMatchAverage();			  
			  
 return true;
}

bool CalculateEquivalentCircularTumor( void )
{
 extern Matrix X;
 extern Matrix Y;
 extern Matrix Phi;
 extern Matrix PhiNecrotic;
 extern Matrix PhiQuiescent;
 extern Matrix Oxygen;
 extern Matrix Mobility;
 
 extern double dx;
 extern double dy;
 
 extern double _EPS;
 extern double _PI;
 
 double TumorArea = 0.0;
 double ProliferatingArea = 0.0;
 double HypoxicArea = 0.0;
 double NecroticArea = 0.0;
 double ExteriorArea =0.0;
 
 double MeanExteriorOxygen = 0.0;
 double MeanInteriorMobility = 0.0;
 
 int ExteriorCount = 0;
 int InteriorCount = 0;
 
 extern double QuiescenceOxygenLevel;
 extern double NecrosisOxygenLevel;
 
 int ProliferatingCount = 0;
 int HypoxicCount = 0;
 int NecroticCount = 0;
 
 // calculate the areas and any appropriate averages, as well as the centroid
 
 for( int i=0 ; i < Phi.TellRows() ; i++ )
 {
  for( int j=0 ; j < Phi.TellCols() ; j++ )
  {
   if( *Phi(i,j) <= _EPS )
   {
    TumorCenterX += *X(i);
	TumorCenterY += *Y(j);
	InteriorCount++;
	
	bool done = false;
	if( *PhiNecrotic(i,j) <= _EPS )
	{
	 NecroticCount++;
	 done = true;
	}
	if( *Oxygen(i,j) < NecrosisOxygenLevel && !done )
	{ 
	 NecroticCount++;
	 done = true; 
	}
	
    if( *PhiQuiescent(i,j) <= _EPS )
    {
     HypoxicCount++;
     done = true;
    }
	if( *Oxygen(i,j) < QuiescenceOxygenLevel && !done )
	{ 
	 HypoxicCount++;
	 done = true; 
	}
	
	if( !done )
	{ ProliferatingCount++; }
 
	MeanInteriorMobility += *Mobility(i,j);
   }
   else
   {
    MeanExteriorOxygen += *Oxygen(i,j);
    ExteriorCount++;
   }
  
  }
 }

 TumorCenterX /= (double) InteriorCount;
 TumorCenterY /= (double) InteriorCount;
 
 MeanInteriorMobility /= (double) InteriorCount;
 mu_in= MeanInteriorMobility;
 mu_out = 1.0;
 
 MeanExteriorOxygen /= (double) ExteriorCount;
 o2_ext = MeanExteriorOxygen;
 o2_H = QuiescenceOxygenLevel;
 o2_N = NecrosisOxygenLevel;
 
 ProliferatingArea = ProliferatingCount*dx*dy;
 HypoxicArea = HypoxicCount*dx*dy;
 NecroticArea = NecroticCount*dx*dy;
 ExteriorArea = ExteriorCount*dx*dy;
 TumorArea = ProliferatingArea + NecroticArea + HypoxicArea; 
 
 extern double ApoptosisRate;        // A in all papers
 extern double TumorAggressivenessG; // G in all papers
 extern double TumorNecrosisRateGN;  // GN in all papers 
 double A = ApoptosisRate;
 double G = TumorAggressivenessG;
 double G_N = TumorNecrosisRateGN;
 
/* 
 cout << "stats: proliferating area: " << ProliferatingArea << endl
      << "       hypoxic area      : " << HypoxicArea << endl
	  << "       necrotic area     : " << NecroticArea << endl
	  << "       total-------------: " << TumorArea << " vs " << InteriorCount*dx*dy << endl
	  << "       exterior area     : " << ExteriorArea << endl
	  << " grand total-------------: " << TumorArea + ExteriorArea << " vs " 
	  << dx*dy*Phi.TellRows()*Phi.TellCols() << endl << endl;
	  
 cout << "exterior o2: " << o2_ext << endl
      << "interior mu: " << mu_in << endl << endl;
	   
 cout << "centroid: ( " << TumorCenterX << " , " << TumorCenterY << " )" << endl
      << "          in domain [" << *X(0) << " , " << *X(X.TellCols()-1) << "]^2" << endl << endl;
*/	  

 // calculate the minimum outer radius to fully contain the computational domain

 double minimum_far_radius = square( *X(0) - TumorCenterX ) + square( *Y(0) - TumorCenterY );
 
 double temp =  square( *X( X.TellCols()-1 ) - TumorCenterX ) + square( *Y(0) - TumorCenterY );
 if( temp > minimum_far_radius )
 { minimum_far_radius = temp; }
 
 temp =  square( *X(X.TellCols()-1) - TumorCenterX ) + square( *Y(Y.TellCols()-1) - TumorCenterY );
 if( temp > minimum_far_radius )
 { minimum_far_radius = temp; }
 
 temp =  square( *X(0) - TumorCenterX ) + square( *Y(Y.TellCols()-1) - TumorCenterY );
 if( temp > minimum_far_radius )
 { minimum_far_radius = temp; }
 
 minimum_far_radius = sqrt( minimum_far_radius );
 
 // calculate equivalent radii;
 
 R = sqrt( TumorArea / _PI );
 R_H = sqrt( (HypoxicArea + NecroticArea) / _PI );
 R_N = sqrt( NecroticArea / _PI );
 R_far = R + FarFieldRadialPressure_DissipationDistance;
 if( R_far < minimum_far_radius )
 { R_far = minimum_far_radius; }
 
 // if no hypoxic or necrotic zones, assume so small that external boundary is 
 // not affected.
 
 if( HypoxicCount == 0 || NecroticCount == 0 )
 {
  ExternalRadialPressure_const = 0.0;
  ExternalRadialPressure_logterm = 0.0;
  
  return false;
 }
 
// cout << "radii: " << R_N << " " << R_H << " " << R << " " << R_far << endl << endl;

 // construct the "matching linear system"
 
 Matrix Pressure_mat = zeros(7,7);
 Matrix Pressure_rhs = zeros(7,7);
 
 double sigma_0 = o2_H - R_H*(o2_ext-o2_H)/(R-R_H);
 double sigma_1 = (o2_ext-o2_H)/(R-R_H);
 
 *Pressure_mat(0,0) = 1.0;
 *Pressure_mat(0,1) = log(R_far);
 *Pressure_rhs(0) = 0;

 *Pressure_mat(1,0) = 0;
 *Pressure_mat(1,1) = mu_out/R;
 *Pressure_mat(1,2) = 0;
 *Pressure_mat(1,3) = -mu_in/R;
 *Pressure_rhs(1) = - (sigma_0-A)*R/2.0 - sigma_1*square(R)/3.0;

 *Pressure_mat(2,0) = -1;
 *Pressure_mat(2,1) = -log(R);
 *Pressure_mat(2,2) = 1;
 *Pressure_mat(2,3) = log(R);
 *Pressure_rhs(2) = 1.0/(G*R) + (  (sigma_0-A)*square(R)/4.0 +  (sigma_1/9)*IntPow(R,3) )/mu_in;

 *Pressure_mat(3,0) = 0;
 *Pressure_mat(3,1) = 0;
 *Pressure_mat(3,2) = 1;
 *Pressure_mat(3,3) = log(R_H);
 *Pressure_mat(3,4) = -1;
 *Pressure_mat(3,5) = -log(R_H);
 *Pressure_rhs(3) = ( (sigma_0-A)*square(R_H)/4.0 + sigma_1*IntPow(R_H,3)/9.0  )/mu_in;

 *Pressure_mat(4,0) = 0;
 *Pressure_mat(4,1) = 0;
 *Pressure_mat(4,2) = 0;
 *Pressure_mat(4,3) = 1.0/R_H;
 *Pressure_mat(4,4) = 0;
 *Pressure_mat(4,5) = -1.0/R_H;
 *Pressure_rhs(4) = ( (sigma_0-A)*R_H/2.0 + sigma_1*square(R_H)/3.0 )/mu_in;

 *Pressure_mat(5,0) = 0;
 *Pressure_mat(5,1) = 0;
 *Pressure_mat(5,2) = 0;
 *Pressure_mat(5,3) = 0;
 *Pressure_mat(5,4) = 1;
 *Pressure_mat(5,5) = log(R_N);
 *Pressure_mat(5,6) = -1;
 *Pressure_rhs(5) = G_N*square(R_N)/(4.0*mu_in);

 *Pressure_mat(6,0) = 0.0;
 *Pressure_mat(6,1) = 0.0;
 *Pressure_mat(6,2) = 0.0;
 *Pressure_mat(6,3) = 0.0;
 *Pressure_mat(6,4) = 0.0;
 *Pressure_mat(6,5) = 1;
 *Pressure_mat(6,6) = 0;
 *Pressure_rhs(6) = G_N*R_N/mu_in;

 Matrix Pressure_coeff = GaussianEliminationPartialPivoting( Pressure_mat, Pressure_rhs  );

 // the solution is now:
 // coeff(0) + coeff(1)*log(r) // exterior
 // coeff(2) + coeff(3)*log(r) - ( (sigma_0-A)/4 * r^2 + sigma_1/9 * r^3 )/mu_in // proliferating
 // coeff(4) + coeff(5)*log(r) // hypoxic
 // coeff(6) + GN/(4*mu_in) * r^2 // necrotic
 
 // using matching conditions, etc. 
 
 ExternalRadialPressure_const = *Pressure_coeff(0);
 ExternalRadialPressure_logterm = *Pressure_coeff(1);
 
 fstream of( "radial_pressure.txt" , ios::out);
  of << "Radial tumor Solution:" << endl 
     << "center: (" << TumorCenterX << "," << TumorCenterY << ")" << endl
	 << "necrotic radius R_N: " << R_N << endl
	 << "hypoxic radius  R_H: " << R_H << endl
	 << "tumor radius      R: " << R << endl
	 << "far-field radius R_far: " << R_far << endl << endl;
	 
 of << "Solution in 0   < r < R_N  : " << endl
    << *Pressure_coeff(6) << " + " << (G_N*0.25/mu_in) << " r^2" << endl
    << "Solution in R_N < r < R_H  : " << endl	
	<< *Pressure_coeff(4) << " + " << *Pressure_coeff(5) << " log(r)" << endl
    << "Solution in R_H < r < R    : " << endl	
	<< *Pressure_coeff(2) << " + " << *Pressure_coeff(3) << " log(r) - "
	<< 1.0/mu_in << " ( " << (sigma_0-A)*0.25 << " r^2 + " << sigma_1/9.0 << " r^3 )" << endl
	<< "Solution in R   < r < R_far: " << endl
	<< *Pressure_coeff(0) << " + " << *Pressure_coeff(1) << " log(r)" << endl << endl;
 
 of.close();
 return true;
}

double FarFieldRadialSolution( double x, double y )
{
 double r = sqrt( square( x - TumorCenterX ) + square( y - TumorCenterY ) );
 return ExternalRadialPressure_const + ExternalRadialPressure_logterm*log(r);
}

double FarFieldRadialSolution_left( double y )
{
 extern Matrix X;
 double x = *X(0);
 return FarFieldRadialSolution(x,y);
}

double FarFieldRadialSolution_right( double y )
{
 extern Matrix X;
 double x = *X(X.TellCols()-1);
 return FarFieldRadialSolution(x,y);
}

double FarFieldRadialSolution_bottom( double x )
{
 extern Matrix Y;
 double y = *Y(0);
 return FarFieldRadialSolution(x,y);
}

double FarFieldRadialSolution_top( double x )
{
 extern Matrix Y;
 double y = *Y(Y.TellCols()-1);
 return FarFieldRadialSolution(x,y);
}


bool UpdateOxyPressure( void )
{
 extern double _EPS;
 
 extern int Xnodes;
 extern int Ynodes;
 extern Matrix X;
 extern Matrix Y;
 extern double dx, dy;
 
 extern Matrix Phi;
 extern Matrix PhiNecrotic;
 extern Matrix PhiQuiescent;
 extern Matrix OxygenDiffusivity;
 extern Matrix BloodVesselIndicator;
 extern Matrix Oxygen;
 extern Matrix Mobility;
 extern Matrix Pressure;
 
 extern double ApoptosisRate;
 extern double TumorNecrosisRateGN;
 extern double TumorProliferationConstant;
 extern char* StringBuffer;
 extern char* LogFile;
 extern Matrix TaxisJump;
 extern bool TopBloodVessel;
 extern bool IsFirstOxygenUpdate;
 
  UpdateOxygenDiffusivity();//HERE

 Matrix PressureSource = zeros(Xnodes,Ynodes);
 Matrix PressureSourceU = zeros(Xnodes,Ynodes);
 
 extern Matrix ECMtaxisCoefficient, ECM;
 
 UpdateECMtaxisCoefficient();
 
 TaxisJump = zeros(Xnodes,Ynodes);
 
 extern double FunctionSmoothingLevel;
 Matrix ECMsmoothed = GaussianFilter(X,Y,ECM,FunctionSmoothingLevel*dx); 
 
  
 UpdateMobility();
 
 
  if( Oxygen.TellRows() != Xnodes || Oxygen.TellCols() != Ynodes )
 { 
  Oxygen = zeros(Xnodes,Ynodes); 
  StartStringBuffer( "Warning: Oxygen previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 if( OxygenDiffusivity.TellRows() != Xnodes || OxygenDiffusivity.TellCols() != Ynodes )
 { 
  OxygenDiffusivity = ones(Xnodes,Ynodes);  
  StartStringBuffer( "Warning: Oxygen diffusivity previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 if( Oxygen.TellRows() != Xnodes || Oxygen.TellCols() != Ynodes )
 {
  Oxygen = zeros(Xnodes,Ynodes);
  StartStringBuffer( "Warning: Oxygen diffusivity previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 if( PhiNecrotic.TellRows() != Xnodes || PhiNecrotic.TellCols() != Ynodes )
 {
  PhiNecrotic = ones(Xnodes,Ynodes);
  StartStringBuffer( "Warning: PhiN diffusivity previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 if( PhiQuiescent.TellRows() != Xnodes || PhiQuiescent.TellCols() != Ynodes )
 {
  PhiQuiescent = ones(Xnodes,Ynodes); 
  StartStringBuffer( "Warning: PhiQ diffusivity previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 if( BloodVesselIndicator.TellRows() != Xnodes || BloodVesselIndicator.TellCols() != Ynodes )
 {
  BloodVesselIndicator = zeros(Xnodes,Ynodes);
  StartStringBuffer( "Warning: Blood Vessel Indicator previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }

 // test line
 PhiQuiescent = ones(Xnodes,Ynodes);

 // end test line

char TopO2conditionType = 'D';
 double (*TopO2BC)(double) = O2_OneFunction;
 if( !TopBloodVessel )
 {
  TopO2conditionType = 'N'; 
  TopO2BC = O2_ZeroFunction;
 }
 
 extern bool BottomBloodVessel;
 char BottomO2conditionType = 'D';
 double (*BottomO2BC)(double) = O2_OneFunction;
 if( !BottomBloodVessel )
 {
  BottomO2conditionType = 'N'; 
  BottomO2BC = O2_ZeroFunction;
 }


 int MaxIterations = 20000;
 double OTolerance = IntPow(dx,4);
 double PTolerance = IntPow(dx,4);
 char bctype = 'D'; // 'N';
 
// bctype = PressureBoundaryConditionType;

 double (*PressureBC_left)(double) = PressureBC;
 double (*PressureBC_right)(double) = PressureBC;
 double (*PressureBC_bottom)(double) = PressureBC;
 double (*PressureBC_top)(double) = PressureBC;

 
 if( UseFarFieldRadialDirichletCondition )
 {
  cout << "Using radial far-field tumor pressure for Dirichlet boundary condition ... " << endl;
  CalculateEquivalentCircularTumor();
  bctype = 'D';
  
  PressureBC_left = FarFieldRadialSolution_left;
  PressureBC_right = FarFieldRadialSolution_right;
  PressureBC_bottom = FarFieldRadialSolution_bottom;
  PressureBC_top = FarFieldRadialSolution_top;
 }
 extern double DRG_Effect;
 extern Matrix TracerD;
 extern double TracerMax;
		  
for( int i=1 ; i < Xnodes-1 ; i++ )
 {
  for( int j=1 ; j < Ynodes-1 ; j++ )
  {
   bool Done = false;
   if( *PhiNecrotic(i,j) <= _EPS )
   {
    *PressureSource(i,j) = -TumorNecrosisRateGN; 
    Done = true; 
   }
   if( *PhiQuiescent(i,j) <= _EPS && !Done )
   { Done = true; }
   if( *Phi(i,j) <= _EPS && !Done )
   {
    *PressureSource(i,j) =  TumorProliferationConstant*( *Oxygen(i,j)*(1-DRG_Effect* *TracerD(i,j)*TracerMax) - ApoptosisRate);
   }
   
   double ECMtaxisXr   = 0.5*( *ECMtaxisCoefficient(i+1,j) + *ECMtaxisCoefficient(i,j) );
   double ECMtaxisXell = 0.5*( *ECMtaxisCoefficient(i-1,j) + *ECMtaxisCoefficient(i,j) );
   double ECMtaxisYr   = 0.5*( *ECMtaxisCoefficient(i,j+1) + *ECMtaxisCoefficient(i,j) );
   double ECMtaxisYell = 0.5*( *ECMtaxisCoefficient(i,j-1) + *ECMtaxisCoefficient(i,j) );
   
   double ECM_taxis_contribution = ( ECMtaxisXr*(*ECMsmoothed(i+1,j)-*ECMsmoothed(i,j)) -
                                     ECMtaxisXell*(*ECMsmoothed(i,j)-*ECMsmoothed(i-1,j))  )/square(dx) 
                                 + ( ECMtaxisYr*(*ECMsmoothed(i,j+1)-*ECMsmoothed(i,j)) - 
								     ECMtaxisYell*(*ECMsmoothed(i,j)-*ECMsmoothed(i,j-1))  )/square(dy);
 								 
   *PressureSourceU(i,j) -= ECM_taxis_contribution; 
   
   GeometricDataPacket GDP = AdaptiveGeometricData(X,Y,Phi,*X(i),*Y(j));

   *TaxisJump(i,j) = *ECMtaxisCoefficient(i,j)* (  GDP.NormalX*((*ECMsmoothed(i+1,j)-*ECMsmoothed(i-1,j))/(2.0*dx))
											      +GDP.NormalY*((*ECMsmoothed(i,j+1)-*ECMsmoothed(i,j-1))/(2.0*dy)) );
  }
 }


   WithOxygenGhostFluidMethod_Generic_MatrixForm(X, Y, Phi, Mobility,
			   PressureSource, PressureSourceU, PressureBC_left,bctype, PressureBC_right,bctype, 
			   PressureBC_bottom,bctype, PressureBC_top,bctype, PressureJumpFunction, PressureNormalJumpFunction,
			   PTolerance, O2_ZeroFunction, 'N', O2_ZeroFunction, 'N', BottomO2BC, BottomO2conditionType, TopO2BC, TopO2conditionType, 
			   OxygenDiffusionConstantFunction, OxygenSourceTimesO2, OxygenSourceRemainder, MaxIterations, OTolerance, ImplicitUpdate, UpdateBoundaryPoint);

			  
 //ShiftPressureToMatchAverage();	
 
 UpdateNecrosisLevelset();

 UpdateQuiescentLevelset();	

/******

	 double v = (pow(5, 0.5) - 0.1)/(pow(5,0.5) - 1);
 
 
  for(int i=0; i< Xnodes; i++)
 {
 for(int j=0; j< Ynodes; j++)
 {
 double r = pow( pow((j-120)*1.0/200,2) + pow((i-100)*1.0/200, 2) , 0.5);
 if(r<= 0.1)
 {
 *Pressure(i, j) = 1;
 }
 else
 {
 *Pressure(i, j) =(v - r)*(v - r)*1.0/pow((v-0.1), 2);

 }
 // *TAF(i, j) = exp(-pow((j)*1.0/200,2)/0.45); 
 //cout<<*TAF(i, j) <<endl;

 }

 
 }





*****************/ 
 return true;
}

bool UpdateOxyWTPressure( void )
{
	extern double _EPS;
	
	extern int Xnodes;
	extern int Ynodes;
	extern Matrix X;
	extern Matrix Y;
	extern double dx, dy;
	
	extern Matrix Phi;
	extern Matrix PhiNecrotic;
	extern Matrix PhiQuiescent;
	extern Matrix OxygenDiffusivity;
	extern Matrix BloodVesselIndicator;
	extern Matrix Oxygen;
	extern Matrix Mobility;
	extern Matrix Pressure;
	extern Matrix Leak;
	extern Matrix Drainage;
	extern Matrix TotalPressure;
	extern Matrix OutPressure;
	extern double ApoptosisRate;
	extern double TumorNecrosisRateGN;
	extern double TumorProliferationConstant;
	extern char* StringBuffer;
	extern char* LogFile;
	extern double ChMobility;
	extern double ChLength;
	extern double MitosisRate;
	extern double VFCell;
	extern double VFWater;
	extern double Lym;
	extern double LymPressure;
	extern Matrix Original_ECM;
	extern bool TopBloodVessel;
	extern bool IsFirstOxygenUpdate;
	extern Matrix HomPressure;
	extern Matrix InterPressure, Permeability, Conduct;
	extern Matrix CollapseClock;
    extern Matrix TaxisJump;
	UpdateOxygenDiffusivity();//HERE
	
	Matrix PressureSource = zeros(Xnodes,Ynodes);
	Matrix PressureSourceU = zeros(Xnodes,Ynodes);
	
	extern Matrix ECMtaxisCoefficient, ECM;
	
	UpdateECMtaxisCoefficient();
	
	TaxisJump = zeros(Xnodes,Ynodes);
	Leak = zeros(Xnodes,Ynodes);
	extern double FunctionSmoothingLevel;
	Matrix ECMsmoothed = GaussianFilter(X,Y,ECM,FunctionSmoothingLevel*dx); 
	
	
	UpdateMobility();
	
	
	if( Oxygen.TellRows() != Xnodes || Oxygen.TellCols() != Ynodes )
	{ 
		Oxygen = zeros(Xnodes,Ynodes); 
		StartStringBuffer( "Warning: Oxygen previously uninitialized. Building ... \n" );
		cout << StringBuffer << endl;
		MakeLogEntry( LogFile, StringBuffer );
	}
	if( OxygenDiffusivity.TellRows() != Xnodes || OxygenDiffusivity.TellCols() != Ynodes )
	{ 
		OxygenDiffusivity = ones(Xnodes,Ynodes);  
		StartStringBuffer( "Warning: Oxygen diffusivity previously uninitialized. Building ... \n" );
		cout << StringBuffer << endl;
		MakeLogEntry( LogFile, StringBuffer );
	}
	if( Oxygen.TellRows() != Xnodes || Oxygen.TellCols() != Ynodes )
	{
		Oxygen = zeros(Xnodes,Ynodes);
		StartStringBuffer( "Warning: Oxygen diffusivity previously uninitialized. Building ... \n" );
		cout << StringBuffer << endl;
		MakeLogEntry( LogFile, StringBuffer );
	}
	if( PhiNecrotic.TellRows() != Xnodes || PhiNecrotic.TellCols() != Ynodes )
	{
		PhiNecrotic = ones(Xnodes,Ynodes);
		StartStringBuffer( "Warning: PhiN diffusivity previously uninitialized. Building ... \n" );
		cout << StringBuffer << endl;
		MakeLogEntry( LogFile, StringBuffer );
	}
	if( PhiQuiescent.TellRows() != Xnodes || PhiQuiescent.TellCols() != Ynodes )
	{
		PhiQuiescent = ones(Xnodes,Ynodes); 
		StartStringBuffer( "Warning: PhiQ diffusivity previously uninitialized. Building ... \n" );
		cout << StringBuffer << endl;
		MakeLogEntry( LogFile, StringBuffer );
	}
	if( BloodVesselIndicator.TellRows() != Xnodes || BloodVesselIndicator.TellCols() != Ynodes )
	{
		BloodVesselIndicator = zeros(Xnodes,Ynodes);
		StartStringBuffer( "Warning: Blood Vessel Indicator previously uninitialized. Building ... \n" );
		cout << StringBuffer << endl;
		MakeLogEntry( LogFile, StringBuffer );
	}
	
	// test line
	PhiQuiescent = ones(Xnodes,Ynodes);
	
	// end test line
	
	char TopO2conditionType = 'D';
	double (*TopO2BC)(double) = O2_OneFunction;
	if( !TopBloodVessel )
	{
		TopO2conditionType = 'N'; 
		TopO2BC = O2_ZeroFunction;
	}
	
	extern bool BottomBloodVessel;
	char BottomO2conditionType = 'D';
	double (*BottomO2BC)(double) = O2_OneFunction;
	if( !BottomBloodVessel )
	{
		BottomO2conditionType = 'N'; 
		BottomO2BC = O2_ZeroFunction;
	}
	
	
	int MaxIterations = 20000;
	double OTolerance = IntPow(dx,3);//4
	double PTolerance = IntPow(dx,3);//4
	double WTolerance = IntPow(dx,2);//3
	char bctype = 'D'; // 'N';
	
	// bctype = PressureBoundaryConditionType;
	
	double (*PressureBC_left)(double) = PressureBC;
	double (*PressureBC_right)(double) = PressureBC;
	double (*PressureBC_bottom)(double) = PressureBC;
	double (*PressureBC_top)(double) = PressureBC;
	
	
	if( UseFarFieldRadialDirichletCondition )
	{
		cout << "Using radial far-field tumor pressure for Dirichlet boundary condition ... " << endl;
		CalculateEquivalentCircularTumor();
		bctype = 'D';
		
		PressureBC_left = FarFieldRadialSolution_left;
		PressureBC_right = FarFieldRadialSolution_right;
		PressureBC_bottom = FarFieldRadialSolution_bottom;
		PressureBC_top = FarFieldRadialSolution_top;
	}
	extern double DRG_Effect;
	extern Matrix TracerD;
	extern double TracerMax;
	
	for( int i=1 ; i < Xnodes-1 ; i++ )
	{
		for( int j=1 ; j < Ynodes-1 ; j++ )
		{
			bool Done = false;
			if( *PhiNecrotic(i,j) <= _EPS )
			{
				*PressureSource(i,j) = -TumorNecrosisRateGN; 
				Done = true; 
			}
			if( *PhiQuiescent(i,j) <= _EPS && !Done )
			{ Done = true; }
			if( *Phi(i,j) <= _EPS && !Done )
			{
				*PressureSource(i,j) =  TumorProliferationConstant*( *Oxygen(i,j)*(1-DRG_Effect* *TracerD(i,j)*TracerMax) - ApoptosisRate);
			}
			
			double ECMtaxisXr   = 0.5*( *ECMtaxisCoefficient(i+1,j) + *ECMtaxisCoefficient(i,j) );
			double ECMtaxisXell = 0.5*( *ECMtaxisCoefficient(i-1,j) + *ECMtaxisCoefficient(i,j) );
			double ECMtaxisYr   = 0.5*( *ECMtaxisCoefficient(i,j+1) + *ECMtaxisCoefficient(i,j) );
			double ECMtaxisYell = 0.5*( *ECMtaxisCoefficient(i,j-1) + *ECMtaxisCoefficient(i,j) );
			
			double ECM_taxis_contribution = ( ECMtaxisXr*(*ECMsmoothed(i+1,j)-*ECMsmoothed(i,j)) -
											 ECMtaxisXell*(*ECMsmoothed(i,j)-*ECMsmoothed(i-1,j))  )/square(dx) 
			+ ( ECMtaxisYr*(*ECMsmoothed(i,j+1)-*ECMsmoothed(i,j)) - 
			   ECMtaxisYell*(*ECMsmoothed(i,j)-*ECMsmoothed(i,j-1))  )/square(dy);
			
			*PressureSourceU(i,j) -= ECM_taxis_contribution; 
			
			GeometricDataPacket GDP = AdaptiveGeometricData(X,Y,Phi,*X(i),*Y(j));
			
			*TaxisJump(i,j) = *ECMtaxisCoefficient(i,j)* (  GDP.NormalX*((*ECMsmoothed(i+1,j)-*ECMsmoothed(i-1,j))/(2.0*dx))
														  +GDP.NormalY*((*ECMsmoothed(i,j+1)-*ECMsmoothed(i,j-1))/(2.0*dy)) );
		}
	}
	
	
	WithOxygenWTGhostFluidMethod_Generic_MatrixForm(X, Y, Phi, Mobility,
												  PressureSource, PressureSourceU, PressureBC_left,bctype, PressureBC_right,bctype, 
												  PressureBC_bottom,bctype, PressureBC_top,bctype, PressureJumpFunction, PressureNormalJumpFunction,
												  PTolerance, O2_ZeroFunction, 'N', O2_ZeroFunction, 'N', BottomO2BC, BottomO2conditionType, TopO2BC, TopO2conditionType, 
													OxygenDiffusionConstantFunction, OxygenWTSourceTimesO2, OxygenWTSourceRemainder, MaxIterations, OTolerance, ImplicitUpdate, UpdateBoundaryPoint,WTolerance,
													 AngioFluidPressureUpdate,
													 FluidPressureUpdateBoundaryPoint);
	
	
	//ShiftPressureToMatchAverage();	
	
	UpdateNecrosisLevelset();
	
	UpdateQuiescentLevelset();	
	OutPressure = zeros(Xnodes,Ynodes);
	Drainage = zeros(Xnodes,Ynodes);
	TotalPressure = zeros(Xnodes,Ynodes);
    double temp=0;
	double vtemp=0;
	extern double waterpressureA;
	extern double vesselpressureA;
	extern double hompressureA;
	extern Matrix BloodNodePressure;
	waterpressureA = 0;
	vesselpressureA = 0;
	hompressureA = 0;
	for(int i=0; i< Xnodes; i++)
	{
		for(int j=0; j< Ynodes; j++)
		{
			
			*OutPressure(i,j) = *Pressure(i,j) * MitosisRate * ChLength* ChLength / ChMobility;
			
			//cout<<*OutPressure(i,j)<<"  ";
			if(*BloodVesselIndicator(i,j)>0)
			{
				//*Leak(i,j) = 1000* *Conduct(i,j)* *Permeability(i,j)*(*HomPressure(i,j)-VFWater* *InterPressure(i,j)-VFCell * *OutPressure(i,j));
				*Leak(i,j) = 1000* *Conduct(i,j)* *Permeability(i,j)*(*HomPressure(i,j)-VFWater* *InterPressure(i,j));
				if(*Leak(i,j)<0)
				{
					*Leak(i,j)= 0;	
				}
				if(*CollapseClock(i,j)>0)
				{
				   *Leak(i,j)= 0;	
				}
				if((*Phi(i,j)<0) && (*BloodNodePressure(i,j)>0))
				{
				vtemp++;
				vesselpressureA += *BloodNodePressure(i,j);
				hompressureA += *HomPressure(i,j);
				}
				
			}
			//*Drainage(i,j) = Lym*(VFWater*(LymPressure-*InterPressure(i,j))-VFCell * *OutPressure(i,j));
			*Drainage(i,j) = - 1000* *Original_ECM(i,j) *Lym*(LymPressure-VFWater* *InterPressure(i,j)-VFCell * *OutPressure(i,j));
			if(*Drainage(i,j)<0)
			{
				*Drainage(i,j) = 0;	
			}
			*TotalPressure(i,j) = VFCell* *OutPressure(i,j) + VFWater * *InterPressure(i,j);
		    if(*Phi(i,j)<0)
			{
			temp++;
			waterpressureA += *InterPressure(i,j);				
							
			}
		
		
		}
		
		
	}
	waterpressureA /=temp;
	vesselpressureA /=vtemp;
	hompressureA /=temp;
	
	
	/******
	 
	 double v = (pow(5, 0.5) - 0.1)/(pow(5,0.5) - 1);
	 
	 
	 for(int i=0; i< Xnodes; i++)
	 {
	 for(int j=0; j< Ynodes; j++)
	 {
	 double r = pow( pow((j-120)*1.0/200,2) + pow((i-100)*1.0/200, 2) , 0.5);
	 if(r<= 0.1)
	 {
	 *Pressure(i, j) = 1;
	 }
	 else
	 {
	 *Pressure(i, j) =(v - r)*(v - r)*1.0/pow((v-0.1), 2);
	 
	 }
	 // *TAF(i, j) = exp(-pow((j)*1.0/200,2)/0.45); 
	 //cout<<*TAF(i, j) <<endl;
	 
	 }
	 
	 
	 }
	 
	 
	 
	 
	 
	 *****************/ 
	return true;
}


